Introduction
Column
Marginal Effects
After running a regression model of some sort, the common way of interpreting the relationships between predictors and the outcome is by interpreting the regression coefficients. In many situations these interpretations are straightforward, however, in some settings the interpretations can be tricky, or the coefficients can simply not be interpreted.
These complications arise most commonly when a model involves an interaction or moderation. In these sort of models, interpretation of the coefficients requires special care due to the relationship between the various interacted predictors, and the set of information obtainable from the coefficients is often limited without resorting to significant algebra.
Marginal effects offer an improvement over simple coefficient interpretation. Marginal means allow you to estimate the average response under a set of conditions; e.g. the average response for each racial group, or the average response when age and weight are at certain values.
Marginal slopes estimate the slope between a predictor and the outcome when other variables are at some specific value; e.g. the average slope on age within each gender, or the average slope on age when weight is at certain values.
In either case, in addition to estimating these marginal effects, we can easily test hypotheses of equality between them via pairwise comparisons. For example, is the average response different by ethnicity, or is the slope on age different between genders.
If you are familiar with interpreting regression coefficients, and specifically interactions, you may have some idea of how to start addressing all the above examples. However, to complete answer these research questions would require more than simply a table of regression coefficients. With marginal effects, one additional set of results can address all these questions.
Finally, marginal effect estimation is a step towards creating informative visualizations of relationships such as interaction plots.
Data and Libraries
Code
Stata
The webuse command can load the directly.
. webuse margex
(Artificial data for margins)
. summarize, sep(0)
Variable | Obs Mean Std. Dev. Min Max
-------------+---------------------------------------------------------
y | 3,000 69.73357 21.53986 0 146.3
outcome | 3,000 .1696667 .3754023 0 1
sex | 3,000 .5006667 .5000829 0 1
group | 3,000 1.828 .7732714 1 3
age | 3,000 39.799 11.54174 20 60
distance | 3,000 58.58566 181.3987 .25 748.8927
ycn | 3,000 74.47263 22.01264 0 153.8
yc | 3,000 .2413333 .4279633 0 1
treatment | 3,000 .5006667 .5000829 0 1
agegroup | 3,000 2.263667 .8316913 1 3
arm | 3,000 1.970667 .7899457 1 3
R (emmeans)
m <- haven::read_dta("http://www.stata-press.com/data/r15/margex.dta")
m <- data.frame(m)
m$sex <- as_factor(m$sex)
m$group <- as_factor(m$group)
head(m) y outcome sex group age distance ycn yc treatment agegroup arm
1 102.6 1 female 1 55 11.75 97.3 1 1 3 1
2 77.2 0 female 1 60 30.75 78.6 0 1 3 1
3 80.5 0 female 1 27 14.25 54.5 0 1 1 1
4 82.5 1 female 1 60 29.25 98.4 1 1 3 1
5 71.6 1 female 1 52 19.00 105.3 1 1 3 1
6 83.2 1 female 1 49 2.50 89.7 0 1 3 1A note on syntax: To make clear which function belongs to which library, I’ve used R’s double colon notation to be explicit. You do not need to include these. E.g. library::function() can be run simply as function(). There’s a bit of possible confusion in that the workhorse function emmeans, belongs to the “emmeans” package, so emmeans::emmeans can be run as just emmeans to the right of the ::.
R (ggeffects)
The ggeffects package requires several other packages to work fully. However, it does not follow proper package framework and requires you to ensure you have the following packages installed as well:
effectsemmeansggplot2
If you do not have those packages loaded, ggeffects will load them on demand (rather than at the time of loading ggeffects). If you do not have a required package installed, ggeffects will instead produce a warning, so keep an eye out for those and install packages as needed.
m <- haven::read_dta("http://www.stata-press.com/data/r15/margex.dta")
m <- data.frame(m)
m$sex <- as_factor(m$sex)
m$group <- as_factor(m$group)
head(m) y outcome sex group age distance ycn yc treatment agegroup arm
1 102.6 1 female 1 55 11.75 97.3 1 1 3 1
2 77.2 0 female 1 60 30.75 78.6 0 1 3 1
3 80.5 0 female 1 27 14.25 54.5 0 1 1 1
4 82.5 1 female 1 60 29.25 98.4 1 1 3 1
5 71.6 1 female 1 52 19.00 105.3 1 1 3 1
6 83.2 1 female 1 49 2.50 89.7 0 1 3 1A note on syntax: To make clear which function belongs to which library, I’ve used R’s double colon notation to be explicit. You do not need to include these. E.g. library::function() can be run simply as function(). There’s a bit of possible confusion in that the workhorse function ggeffect, belongs to the “ggeffects” package, so ggeffects::ggeffect can be run as just ggeffect.
Categorical Variables
Code
Stata
. regress y i.group
Source | SS df MS Number of obs = 3,000
-------------+---------------------------------- F(2, 2997) = 15.48
Model | 14229.0349 2 7114.51743 Prob > F = 0.0000
Residual | 1377203.97 2,997 459.527518 R-squared = 0.0102
-------------+---------------------------------- Adj R-squared = 0.0096
Total | 1391433.01 2,999 463.965657 Root MSE = 21.437
------------------------------------------------------------------------------
y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
group |
2 | .4084991 .8912269 0.46 0.647 -1.338979 2.155977
3 | 5.373138 1.027651 5.23 0.000 3.358165 7.38811
|
_cons | 68.35805 .6190791 110.42 0.000 67.14419 69.57191
------------------------------------------------------------------------------
. margins group
Adjusted predictions Number of obs = 3,000
Model VCE : OLS
Expression : Linear prediction, predict()
------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
group |
1 | 68.35805 .6190791 110.42 0.000 67.14419 69.57191
2 | 68.76655 .6411134 107.26 0.000 67.50948 70.02361
3 | 73.73119 .8202484 89.89 0.000 72.12288 75.33949
------------------------------------------------------------------------------
R (emmeans)
Call:
lm(formula = y ~ group, data = m)
Residuals:
Min 1Q Median 3Q Max
-73.531 -14.758 0.101 14.536 72.569
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 68.3580 0.6191 110.419 < 2e-16 ***
group2 0.4085 0.8912 0.458 0.647
group3 5.3731 1.0277 5.229 1.83e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21.44 on 2997 degrees of freedom
Multiple R-squared: 0.01023, Adjusted R-squared: 0.009566
F-statistic: 15.48 on 2 and 2997 DF, p-value: 2.045e-07 group emmean SE df lower.CL upper.CL
1 68.4 0.619 2997 67.1 69.6
2 68.8 0.641 2997 67.5 70.0
3 73.7 0.820 2997 72.1 75.3
Confidence level used: 0.95 R (ggeffects)
Call:
lm(formula = y ~ group, data = m)
Residuals:
Min 1Q Median 3Q Max
-73.531 -14.758 0.101 14.536 72.569
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 68.3580 0.6191 110.419 < 2e-16 ***
group2 0.4085 0.8912 0.458 0.647
group3 5.3731 1.0277 5.229 1.83e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21.44 on 2997 degrees of freedom
Multiple R-squared: 0.01023, Adjusted R-squared: 0.009566
F-statistic: 15.48 on 2 and 2997 DF, p-value: 2.045e-07
# Predicted values of y
# x = group
x predicted std.error conf.low conf.high
1 68.358 0.619 67.144 69.572
2 68.767 0.641 67.509 70.024
3 73.731 0.820 72.123 75.339Note that we could have run that as ggeffects::ggeffect(mod1, ~ group), with a formula. However, when we later want to specify certain values of the predictors, we need to use the character version, so I stick with it here.
Math
Assume a linear regression set up with a single categorical variable, \(G\), with three groups. We fit
\[ E(Y|G) = \beta_0 + \beta_1g_2 + \beta_2g_3, \]
where \(g_2\) and \(g_3\) are dummy variables representing membership in groups 2 and 3 respectively (\(g_{2i} = 1\) if observation \(i\) has \(G = 2\), and equal to 0 otherwise.) Since \(g_1\) is not found in the model, it is the reference category.
Therefore, \(\beta_0\) represents the average response among the reference category \(G = 1\). \(\beta_1\) represents the difference in the average response between groups \(G = 1\) and \(G = 2\). Therefore, \(\beta_0 + \beta_1\) is the average response in group \(G = 2\). A similar argument can be made about \(\beta_2\) and group 3.
Interactions between Categorical Variables
Code
Stata
. regress y i.sex##i.group
Source | SS df MS Number of obs = 3,000
-------------+---------------------------------- F(5, 2994) = 92.91
Model | 186898.199 5 37379.6398 Prob > F = 0.0000
Residual | 1204534.81 2,994 402.316235 R-squared = 0.1343
-------------+---------------------------------- Adj R-squared = 0.1329
Total | 1391433.01 2,999 463.965657 Root MSE = 20.058
------------------------------------------------------------------------------
y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
sex |
female | 21.62507 1.509999 14.32 0.000 18.66433 24.58581
|
group |
2 | 11.37444 1.573314 7.23 0.000 8.289552 14.45932
3 | 21.6188 1.588487 13.61 0.000 18.50416 24.73344
|
sex#group |
female#2 | -4.851582 1.942744 -2.50 0.013 -8.66083 -1.042333
female#3 | -6.084875 3.004638 -2.03 0.043 -11.97624 -.1935115
|
_cons | 50.6107 1.367932 37.00 0.000 47.92852 53.29288
------------------------------------------------------------------------------
. margins sex#group
Adjusted predictions Number of obs = 3,000
Model VCE : OLS
Expression : Linear prediction, predict()
------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
sex#group |
male#1 | 50.6107 1.367932 37.00 0.000 47.92852 53.29288
male#2 | 61.98514 .7772248 79.75 0.000 60.46119 63.50908
male#3 | 72.2295 .8074975 89.45 0.000 70.64619 73.8128
female#1 | 72.23577 .63942 112.97 0.000 70.98203 73.48952
female#2 | 78.75863 .9434406 83.48 0.000 76.90877 80.60849
female#3 | 87.7697 2.468947 35.55 0.000 82.92869 92.6107
------------------------------------------------------------------------------
R (emmeans)
Call:
lm(formula = y ~ sex * group, data = m)
Residuals:
Min 1Q Median 3Q Max
-72.029 -13.285 0.464 13.741 67.964
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 50.611 1.368 36.998 < 2e-16 ***
sexfemale 21.625 1.510 14.321 < 2e-16 ***
group2 11.374 1.573 7.230 6.12e-13 ***
group3 21.619 1.588 13.610 < 2e-16 ***
sexfemale:group2 -4.852 1.943 -2.497 0.0126 *
sexfemale:group3 -6.085 3.005 -2.025 0.0429 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 20.06 on 2994 degrees of freedom
Multiple R-squared: 0.1343, Adjusted R-squared: 0.1329
F-statistic: 92.91 on 5 and 2994 DF, p-value: < 2.2e-16 group sex emmean SE df lower.CL upper.CL
1 male 50.6 1.368 2994 47.9 53.3
2 male 62.0 0.777 2994 60.5 63.5
3 male 72.2 0.807 2994 70.6 73.8
1 female 72.2 0.639 2994 71.0 73.5
2 female 78.8 0.943 2994 76.9 80.6
3 female 87.8 2.469 2994 82.9 92.6
Confidence level used: 0.95 R (ggeffects)
Call:
lm(formula = y ~ sex * group, data = m)
Residuals:
Min 1Q Median 3Q Max
-72.029 -13.285 0.464 13.741 67.964
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 50.611 1.368 36.998 < 2e-16 ***
sexfemale 21.625 1.510 14.321 < 2e-16 ***
group2 11.374 1.573 7.230 6.12e-13 ***
group3 21.619 1.588 13.610 < 2e-16 ***
sexfemale:group2 -4.852 1.943 -2.497 0.0126 *
sexfemale:group3 -6.085 3.005 -2.025 0.0429 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 20.06 on 2994 degrees of freedom
Multiple R-squared: 0.1343, Adjusted R-squared: 0.1329
F-statistic: 92.91 on 5 and 2994 DF, p-value: < 2.2e-16
# Predicted values of y
# x = group
# sex = male
x predicted std.error conf.low conf.high
1 50.611 1.368 47.929 53.293
2 61.985 0.777 60.461 63.509
3 72.229 0.807 70.646 73.813
# sex = female
x predicted std.error conf.low conf.high
1 72.236 0.639 70.982 73.490
2 78.759 0.943 76.909 80.608
3 87.770 2.469 82.929 92.611Math
Consider the simplest case, with two binary variables \(Z\) and \(K\). We fit the model
\[ E(Y|Z,K) = \beta_0 + \beta_1Z + \beta_2K + \beta_3ZK, \]
where \(ZK = 0\) only if both \(Z = 1\) and \(K = 1\).
This is functionally equivalent to defining a new variable \(L\),
\[ L = \begin{cases} 0, & K = 0 \& Z = 0 \\ 1, & K = 0 \& Z = 1 \\ 2, & K = 1 \& Z = 0 \\ 3, & K = 1 \& Z = 1, \end{cases} \]
and fitting the model
\[ E(Y|L) = \beta_0 + \beta_1l_1 + \beta_2l_2 + \beta_3l_3. \]
Continuous Variables
Code
Stata
. regress y i.sex c.age c.distance
Source | SS df MS Number of obs = 3,000
-------------+---------------------------------- F(3, 2996) = 143.27
Model | 174576.269 3 58192.0896 Prob > F = 0.0000
Residual | 1216856.74 2,996 406.16046 R-squared = 0.1255
-------------+---------------------------------- Adj R-squared = 0.1246
Total | 1391433.01 2,999 463.965657 Root MSE = 20.153
------------------------------------------------------------------------------
y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
sex |
female | 13.98753 .7768541 18.01 0.000 12.46431 15.51075
age | -.5038264 .0336538 -14.97 0.000 -.5698133 -.4378394
distance | -.0060725 .0020291 -2.99 0.003 -.0100511 -.0020939
_cons | 83.13802 1.327228 62.64 0.000 80.53565 85.74039
------------------------------------------------------------------------------
. margins, at(age = (30 40))
Predictive margins Number of obs = 3,000
Model VCE : OLS
Expression : Linear prediction, predict()
1._at : age = 30
2._at : age = 40
------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_at |
1 | 74.67056 .4941029 151.12 0.000 73.70175 75.63938
2 | 69.6323 .3680117 189.21 0.000 68.91072 70.35388
------------------------------------------------------------------------------
. margins, at(age = (30 40) distance = (10 100))
Predictive margins Number of obs = 3,000
Model VCE : OLS
Expression : Linear prediction, predict()
1._at : age = 30
distance = 10
2._at : age = 30
distance = 100
3._at : age = 40
distance = 10
4._at : age = 40
distance = 100
------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_at |
1 | 74.9656 .5034959 148.89 0.000 73.97837 75.95283
2 | 74.41907 .5014946 148.39 0.000 73.43576 75.40238
3 | 69.92733 .3809974 183.54 0.000 69.18029 70.67438
4 | 69.38081 .3774763 183.80 0.000 68.64067 70.12095
------------------------------------------------------------------------------
. margins sex, at(age = (30 40))
Predictive margins Number of obs = 3,000
Model VCE : OLS
Expression : Linear prediction, predict()
1._at : age = 30
2._at : age = 40
------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_at#sex |
1#male | 67.66747 .5597763 120.88 0.000 66.56989 68.76506
1#female | 81.655 .6902601 118.30 0.000 80.30157 83.00843
2#male | 62.62921 .5370234 116.62 0.000 61.57624 63.68218
2#female | 76.61674 .5331296 143.71 0.000 75.5714 77.66208
------------------------------------------------------------------------------
R (emmeans)
Note that any variable specified in the at option must also appear in the formula. Any variables you place in the formula but not at will be examined at each unique level (so don’t put a continuous variable there!).
Call:
lm(formula = y ~ sex + age + distance, data = m)
Residuals:
Min 1Q Median 3Q Max
-72.356 -13.792 -0.275 13.624 74.950
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 83.138023 1.327228 62.640 < 2e-16 ***
sexfemale 13.987531 0.776854 18.005 < 2e-16 ***
age -0.503826 0.033654 -14.971 < 2e-16 ***
distance -0.006073 0.002029 -2.993 0.00279 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 20.15 on 2996 degrees of freedom
Multiple R-squared: 0.1255, Adjusted R-squared: 0.1246
F-statistic: 143.3 on 3 and 2996 DF, p-value: < 2.2e-16 age emmean SE df lower.CL upper.CL
30 74.7 0.494 2996 73.7 75.6
40 69.6 0.368 2996 68.9 70.3
Results are averaged over the levels of: sex
Confidence level used: 0.95 age distance emmean SE df lower.CL upper.CL
30 10 75.0 0.503 2996 74.0 75.9
40 10 69.9 0.381 2996 69.2 70.7
30 100 74.4 0.501 2996 73.4 75.4
40 100 69.4 0.377 2996 68.6 70.1
Results are averaged over the levels of: sex
Confidence level used: 0.95 age sex emmean SE df lower.CL upper.CL
30 male 67.7 0.560 2996 66.6 68.8
40 male 62.6 0.537 2996 61.6 63.7
30 female 81.7 0.690 2996 80.3 83.0
40 female 76.6 0.533 2996 75.6 77.7
Confidence level used: 0.95 R (ggeffects)
Note the particular syntax for specifying particular values - a quoted string containing the name of the variable, followed by a list in square brackets of comma-separated values (a single value would be just [4]). Whitespace is ignored.
Call:
lm(formula = y ~ sex + age + distance, data = m)
Residuals:
Min 1Q Median 3Q Max
-72.356 -13.792 -0.275 13.624 74.950
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 83.138023 1.327228 62.640 < 2e-16 ***
sexfemale 13.987531 0.776854 18.005 < 2e-16 ***
age -0.503826 0.033654 -14.971 < 2e-16 ***
distance -0.006073 0.002029 -2.993 0.00279 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 20.15 on 2996 degrees of freedom
Multiple R-squared: 0.1255, Adjusted R-squared: 0.1246
F-statistic: 143.3 on 3 and 2996 DF, p-value: < 2.2e-16
# Predicted values of y
# x = age
x predicted std.error conf.low conf.high
30 74.671 0.494 73.702 75.639
40 69.632 0.368 68.911 70.354
# Predicted values of y
# x = age
# distance = 10
x predicted std.error conf.low conf.high
30 74.966 0.503 73.978 75.953
40 69.927 0.381 69.180 70.674
# distance = 100
x predicted std.error conf.low conf.high
30 74.419 0.501 73.436 75.402
40 69.381 0.377 68.641 70.121
# Predicted values of y
# x = sex
# age = 30
x predicted std.error conf.low conf.high
male 67.667 0.56 66.570 68.765
female 81.655 0.69 80.302 83.008
# age = 40
x predicted std.error conf.low conf.high
male 62.629 0.537 61.576 63.682
female 76.617 0.533 75.571 77.662Math
Consider a model with a single continuous predictor, plus a control variable (which may be continuous or categorical).
\[ E(Y|X, Z) = \beta_0 + \beta_1X + \beta_2Z. \]
For any given value of \(X\), it is easy to compute the marginal mean for a given individual. For example,
\[ E(Y|X = 2, Z) = \beta_0 + 2\beta_1 + \beta_2Z. \]
Therefore we can easily compute \(E(y_i|x = 2, z = z_i)\) for each individual (the predicted outcome if each individual had a value of \(X = 2\), with their other values at the observed value) and average them.
Categorical and Continuous Interactions
Code
Stata
. regress y i.sex##c.age
Source | SS df MS Number of obs = 3,000
-------------+---------------------------------- F(3, 2996) = 139.91
Model | 170983.675 3 56994.5583 Prob > F = 0.0000
Residual | 1220449.33 2,996 407.35959 R-squared = 0.1229
-------------+---------------------------------- Adj R-squared = 0.1220
Total | 1391433.01 2,999 463.965657 Root MSE = 20.183
------------------------------------------------------------------------------
y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
sex |
female | 14.92308 2.789012 5.35 0.000 9.454508 20.39165
age | -.4929608 .0480944 -10.25 0.000 -.5872622 -.3986595
|
sex#c.age |
female | -.0224116 .0674167 -0.33 0.740 -.1545994 .1097762
|
_cons | 82.36936 1.812958 45.43 0.000 78.8146 85.92413
------------------------------------------------------------------------------
. margins, at(age = (30 40))
Predictive margins Number of obs = 3,000
Model VCE : OLS
Expression : Linear prediction, predict()
1._at : age = 30
2._at : age = 40
------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_at |
1 | 74.71541 .5089298 146.81 0.000 73.71752 75.71329
2 | 69.67359 .3890273 179.10 0.000 68.9108 70.43638
------------------------------------------------------------------------------
. margins sex, at(age = (30 40))
Adjusted predictions Number of obs = 3,000
Model VCE : OLS
Expression : Linear prediction, predict()
1._at : age = 30
2._at : age = 40
------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_at#sex |
1#male | 67.58054 .5984013 112.94 0.000 66.40722 68.75386
1#female | 81.83127 .8228617 99.45 0.000 80.21784 83.4447
2#male | 62.65093 .5541361 113.06 0.000 61.56441 63.73746
2#female | 76.67755 .5461908 140.39 0.000 75.6066 77.74849
------------------------------------------------------------------------------
R (emmeans)
Call:
lm(formula = y ~ sex * age, data = m)
Residuals:
Min 1Q Median 3Q Max
-71.817 -13.848 -0.116 13.758 75.263
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 82.36936 1.81296 45.434 < 2e-16 ***
sexfemale 14.92308 2.78901 5.351 9.42e-08 ***
age -0.49296 0.04809 -10.250 < 2e-16 ***
sexfemale:age -0.02241 0.06742 -0.332 0.74
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 20.18 on 2996 degrees of freedom
Multiple R-squared: 0.1229, Adjusted R-squared: 0.122
F-statistic: 139.9 on 3 and 2996 DF, p-value: < 2.2e-16 age emmean SE df lower.CL upper.CL
30 74.7 0.509 2996 73.7 75.7
40 69.7 0.389 2996 68.9 70.4
Results are averaged over the levels of: sex
Confidence level used: 0.95 age sex emmean SE df lower.CL upper.CL
30 male 67.6 0.598 2996 66.4 68.8
40 male 62.7 0.554 2996 61.6 63.7
30 female 81.8 0.823 2996 80.2 83.4
40 female 76.7 0.546 2996 75.6 77.7
Confidence level used: 0.95 R (ggeffects)
Call:
lm(formula = y ~ sex * age, data = m)
Residuals:
Min 1Q Median 3Q Max
-71.817 -13.848 -0.116 13.758 75.263
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 82.36936 1.81296 45.434 < 2e-16 ***
sexfemale 14.92308 2.78901 5.351 9.42e-08 ***
age -0.49296 0.04809 -10.250 < 2e-16 ***
sexfemale:age -0.02241 0.06742 -0.332 0.74
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 20.18 on 2996 degrees of freedom
Multiple R-squared: 0.1229, Adjusted R-squared: 0.122
F-statistic: 139.9 on 3 and 2996 DF, p-value: < 2.2e-16
# Predicted values of y
# x = age
x predicted std.error conf.low conf.high
30 74.715 0.509 73.718 75.713
40 69.674 0.389 68.911 70.436
# Predicted values of y
# x = sex
# age = 30
x predicted std.error conf.low conf.high
male 67.581 0.598 66.407 68.754
female 81.831 0.823 80.218 83.445
# age = 40
x predicted std.error conf.low conf.high
male 62.651 0.554 61.564 63.737
female 76.678 0.546 75.607 77.748Marginal Slopes
Code
Stata
. regress y i.sex c.age
Source | SS df MS Number of obs = 3,000
-------------+---------------------------------- F(2, 2997) = 209.88
Model | 170938.657 2 85469.3283 Prob > F = 0.0000
Residual | 1220494.35 2,997 407.238688 R-squared = 0.1229
-------------+---------------------------------- Adj R-squared = 0.1223
Total | 1391433.01 2,999 463.965657 Root MSE = 20.18
------------------------------------------------------------------------------
y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
sex |
female | 14.03271 .7777377 18.04 0.000 12.50775 15.55766
age | -.5043666 .033698 -14.97 0.000 -.5704402 -.4382931
_cons | 82.78115 1.323613 62.54 0.000 80.18586 85.37643
------------------------------------------------------------------------------
. margins, dydx(age)
Average marginal effects Number of obs = 3,000
Model VCE : OLS
Expression : Linear prediction, predict()
dy/dx w.r.t. : age
------------------------------------------------------------------------------
| Delta-method
| dy/dx Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
age | -.5043666 .033698 -14.97 0.000 -.5704402 -.4382931
------------------------------------------------------------------------------
R (emmeans)
Call:
lm(formula = y ~ sex + age, data = m)
Residuals:
Min 1Q Median 3Q Max
-71.99 -13.88 -0.15 13.76 75.28
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 82.7811 1.3236 62.54 <2e-16 ***
sexfemale 14.0327 0.7777 18.04 <2e-16 ***
age -0.5044 0.0337 -14.97 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 20.18 on 2997 degrees of freedom
Multiple R-squared: 0.1229, Adjusted R-squared: 0.1223
F-statistic: 209.9 on 2 and 2997 DF, p-value: < 2.2e-16 1 age.trend SE df lower.CL upper.CL
overall -0.504 0.0337 2997 -0.57 -0.438
Results are averaged over the levels of: sex
Confidence level used: 0.95 Math
The Stata margins command to obtain marginal means includes the “dydx” option, which points to derivative - and indeed, this is exactly what is computer. If we have a basic regression model,
\[ E(Y|X) = \beta_0 + \beta_1X \]
taking the derivative with respect to \(X\) will obtain \(\beta_1\), which is the estimated coefficient.
If \(X\) enters the model in a more complex way, say a polynomial term:
\[ E(Y|X) = \beta_0 + \beta_1X + \beta_2X^2 \]
now the derivative is \(\beta_1 + 2\beta_2X\). Similar to marginal means, where the predicted outcome was estimated for each individual then those outcomes averaged, here the derivative is estimated plugging in each observation’s value of \(X\), then averaged.
Categorical by continuous interaction
Code
Stata
. regress y i.sex##c.age
Source | SS df MS Number of obs = 3,000
-------------+---------------------------------- F(3, 2996) = 139.91
Model | 170983.675 3 56994.5583 Prob > F = 0.0000
Residual | 1220449.33 2,996 407.35959 R-squared = 0.1229
-------------+---------------------------------- Adj R-squared = 0.1220
Total | 1391433.01 2,999 463.965657 Root MSE = 20.183
------------------------------------------------------------------------------
y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
sex |
female | 14.92308 2.789012 5.35 0.000 9.454508 20.39165
age | -.4929608 .0480944 -10.25 0.000 -.5872622 -.3986595
|
sex#c.age |
female | -.0224116 .0674167 -0.33 0.740 -.1545994 .1097762
|
_cons | 82.36936 1.812958 45.43 0.000 78.8146 85.92413
------------------------------------------------------------------------------
. margins sex, dydx(age)
Average marginal effects Number of obs = 3,000
Model VCE : OLS
Expression : Linear prediction, predict()
dy/dx w.r.t. : age
------------------------------------------------------------------------------
| Delta-method
| dy/dx Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
age |
sex |
male | -.4929608 .0480944 -10.25 0.000 -.5872622 -.3986595
female | -.5153724 .0472435 -10.91 0.000 -.6080054 -.4227395
------------------------------------------------------------------------------
R (emmeans)
Call:
lm(formula = y ~ sex * age, data = m)
Residuals:
Min 1Q Median 3Q Max
-71.817 -13.848 -0.116 13.758 75.263
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 82.36936 1.81296 45.434 < 2e-16 ***
sexfemale 14.92308 2.78901 5.351 9.42e-08 ***
age -0.49296 0.04809 -10.250 < 2e-16 ***
sexfemale:age -0.02241 0.06742 -0.332 0.74
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 20.18 on 2996 degrees of freedom
Multiple R-squared: 0.1229, Adjusted R-squared: 0.122
F-statistic: 139.9 on 3 and 2996 DF, p-value: < 2.2e-16 sex age.trend SE df lower.CL upper.CL
male -0.493 0.0481 2996 -0.587 -0.399
female -0.515 0.0472 2996 -0.608 -0.423
Confidence level used: 0.95 Math
Let’s assume we have a binary variable, \(Z\), interacted with a continuous variable, \(X\).
\[ E(Y|X,Z) = \beta_0 + \beta_1X + \beta_2Z + \beta_3XZ \]
Here we’re combining the math from marginal effects with marginal slopes. First, we generate two equations, one for \(Z = 0\) and one for \(Z = 1\):
\[\begin{align*} E(Y|X, Z = 0) & = \beta_0 + \beta_1X \\ E(Y|X, Z = 1) & = \beta_0 + \beta_1X + \beta_2 + \beta_3X = (\beta_0 + \beta_2) + (\beta_1 + \beta_3)X. \end{align*}\]
Then when we differentiate each with respect to \(X\), we obtain \(\beta_1\) and \((\beta_1 + \beta_3)\) respectively.
Categorical Variables
Code
Stata
. regress y i.group
Source | SS df MS Number of obs = 3,000
-------------+---------------------------------- F(2, 2997) = 15.48
Model | 14229.0349 2 7114.51743 Prob > F = 0.0000
Residual | 1377203.97 2,997 459.527518 R-squared = 0.0102
-------------+---------------------------------- Adj R-squared = 0.0096
Total | 1391433.01 2,999 463.965657 Root MSE = 21.437
------------------------------------------------------------------------------
y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
group |
2 | .4084991 .8912269 0.46 0.647 -1.338979 2.155977
3 | 5.373138 1.027651 5.23 0.000 3.358165 7.38811
|
_cons | 68.35805 .6190791 110.42 0.000 67.14419 69.57191
------------------------------------------------------------------------------
. margins group, pwcompare(pv)
Pairwise comparisons of adjusted predictions Number of obs = 3,000
Model VCE : OLS
Expression : Linear prediction, predict()
-----------------------------------------------------
| Delta-method Unadjusted
| Contrast Std. Err. t P>|t|
-------------+---------------------------------------
group |
2 vs 1 | .4084991 .8912269 0.46 0.647
3 vs 1 | 5.373138 1.027651 5.23 0.000
3 vs 2 | 4.964638 1.041073 4.77 0.000
-----------------------------------------------------
R (emmeans)
Call:
lm(formula = y ~ group, data = m)
Residuals:
Min 1Q Median 3Q Max
-73.531 -14.758 0.101 14.536 72.569
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 68.3580 0.6191 110.419 < 2e-16 ***
group2 0.4085 0.8912 0.458 0.647
group3 5.3731 1.0277 5.229 1.83e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21.44 on 2997 degrees of freedom
Multiple R-squared: 0.01023, Adjusted R-squared: 0.009566
F-statistic: 15.48 on 2 and 2997 DF, p-value: 2.045e-07Note that pairs is a generic, meaning that emmeans::pairs is an invalid call. Despite this, the “emmeans” package is required to call the below.
contrast estimate SE df t.ratio p.value
1 - 2 -0.408 0.891 2997 -0.458 0.6467
1 - 3 -5.373 1.028 2997 -5.229 <.0001
2 - 3 -4.965 1.041 2997 -4.769 <.0001 The adjust = "none" argument skips any multiple comparisons correction. This is done to obtain the same results as the regression coefficients. If you’d prefer a more ANOVA post-hoc approach, there are other options for adjust, see ?emmeans::contrast for details.
Categorical Interactions
Column
Stata
. regress y i.group##i.sex
Source | SS df MS Number of obs = 3,000
-------------+---------------------------------- F(5, 2994) = 92.91
Model | 186898.199 5 37379.6398 Prob > F = 0.0000
Residual | 1204534.81 2,994 402.316235 R-squared = 0.1343
-------------+---------------------------------- Adj R-squared = 0.1329
Total | 1391433.01 2,999 463.965657 Root MSE = 20.058
------------------------------------------------------------------------------
y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
group |
2 | 11.37444 1.573314 7.23 0.000 8.289552 14.45932
3 | 21.6188 1.588487 13.61 0.000 18.50416 24.73344
|
sex |
female | 21.62507 1.509999 14.32 0.000 18.66433 24.58581
|
group#sex |
2#female | -4.851582 1.942744 -2.50 0.013 -8.66083 -1.042333
3#female | -6.084875 3.004638 -2.03 0.043 -11.97624 -.1935115
|
_cons | 50.6107 1.367932 37.00 0.000 47.92852 53.29288
------------------------------------------------------------------------------
. margins group#sex, pwcompare(pv)
Pairwise comparisons of adjusted predictions Number of obs = 3,000
Model VCE : OLS
Expression : Linear prediction, predict()
------------------------------------------------------------------
| Delta-method Unadjusted
| Contrast Std. Err. t P>|t|
--------------------------+---------------------------------------
group#sex |
(1#female) vs (1#male) | 21.62507 1.509999 14.32 0.000
(2#male) vs (1#male) | 11.37444 1.573314 7.23 0.000
(2#female) vs (1#male) | 28.14793 1.661722 16.94 0.000
(3#male) vs (1#male) | 21.6188 1.588487 13.61 0.000
(3#female) vs (1#male) | 37.159 2.822577 13.16 0.000
(2#male) vs (1#female) | -10.25064 1.006447 -10.18 0.000
(2#female) vs (1#female) | 6.522856 1.13971 5.72 0.000
(3#male) vs (1#female) | -.0062748 1.030005 -0.01 0.995
(3#female) vs (1#female) | 15.53392 2.550404 6.09 0.000
(2#female) vs (2#male) | 16.77349 1.222358 13.72 0.000
(3#male) vs (2#male) | 10.24436 1.120772 9.14 0.000
(3#female) vs (2#male) | 25.78456 2.588393 9.96 0.000
(3#male) vs (2#female) | -6.529131 1.241826 -5.26 0.000
(3#female) vs (2#female) | 9.011069 2.643063 3.41 0.001
(3#female) vs (3#male) | 15.5402 2.597644 5.98 0.000
------------------------------------------------------------------
This produces a lot of comparisons; if we only cared about the comparisons one way (e.g. gender differences within each group), we can restrict our output to these results.
. margins sex@group, contrast(pv nowald)
Contrasts of adjusted predictions Number of obs = 3,000
Model VCE : OLS
Expression : Linear prediction, predict()
------------------------------------------------------------
| Delta-method
| Contrast Std. Err. t P>|t|
--------------------+---------------------------------------
sex@group |
(female vs base) 1 | 21.62507 1.509999 14.32 0.000
(female vs base) 2 | 16.77349 1.222358 13.72 0.000
(female vs base) 3 | 15.5402 2.597644 5.98 0.000
------------------------------------------------------------
Note that the results are identical, it’s just a nicer output.
The pv option requests p-values instead of confidence intervals; the nowald option suppresses an omnibus test of any difference between groups (e.g. ANOVA output).
R (emmeans)
Call:
lm(formula = y ~ group * sex, data = m)
Residuals:
Min 1Q Median 3Q Max
-72.029 -13.285 0.464 13.741 67.964
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 50.611 1.368 36.998 < 2e-16 ***
group2 11.374 1.573 7.230 6.12e-13 ***
group3 21.619 1.588 13.610 < 2e-16 ***
sexfemale 21.625 1.510 14.321 < 2e-16 ***
group2:sexfemale -4.852 1.943 -2.497 0.0126 *
group3:sexfemale -6.085 3.005 -2.025 0.0429 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 20.06 on 2994 degrees of freedom
Multiple R-squared: 0.1343, Adjusted R-squared: 0.1329
F-statistic: 92.91 on 5 and 2994 DF, p-value: < 2.2e-16 contrast estimate SE df t.ratio p.value
male,1 - female,1 -21.62507 1.51 2994 -14.321 <.0001
male,1 - male,2 -11.37444 1.57 2994 -7.230 <.0001
male,1 - female,2 -28.14793 1.66 2994 -16.939 <.0001
male,1 - male,3 -21.61880 1.59 2994 -13.610 <.0001
male,1 - female,3 -37.15900 2.82 2994 -13.165 <.0001
female,1 - male,2 10.25064 1.01 2994 10.185 <.0001
female,1 - female,2 -6.52286 1.14 2994 -5.723 <.0001
female,1 - male,3 0.00627 1.03 2994 0.006 0.9951
female,1 - female,3 -15.53392 2.55 2994 -6.091 <.0001
male,2 - female,2 -16.77349 1.22 2994 -13.722 <.0001
male,2 - male,3 -10.24436 1.12 2994 -9.140 <.0001
male,2 - female,3 -25.78456 2.59 2994 -9.962 <.0001
female,2 - male,3 6.52913 1.24 2994 5.258 <.0001
female,2 - female,3 -9.01107 2.64 2994 -3.409 0.0007
male,3 - female,3 -15.54020 2.60 2994 -5.982 <.0001 The adjust = "none" argument skips any multiple comparisons correction. This is done to obtain the same results as the regression coefficients. If you’d prefer a more ANOVA post-hoc approach, there are other options for adjust, see ?emmeans::contrast for details.
This produces a lot of comparisons; if we only cared about the comparisons one way (e.g. gender differences within each group), we can restrict our output to these results.
group = 1:
contrast estimate SE df t.ratio p.value
male - female -21.6 1.51 2994 -14.321 <.0001
group = 2:
contrast estimate SE df t.ratio p.value
male - female -16.8 1.22 2994 -13.722 <.0001
group = 3:
contrast estimate SE df t.ratio p.value
male - female -15.5 2.60 2994 -5.982 <.0001 Note that the results are identical, it’s just a nicer output.
Marginal Slopes
Column
Stata
. regress y i.group##c.age
Source | SS df MS Number of obs = 3,000
-------------+---------------------------------- F(5, 2994) = 21.15
Model | 47475.1592 5 9495.03183 Prob > F = 0.0000
Residual | 1343957.85 2,994 448.883716 R-squared = 0.0341
-------------+---------------------------------- Adj R-squared = 0.0325
Total | 1391433.01 2,999 463.965657 Root MSE = 21.187
------------------------------------------------------------------------------
y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
group |
2 | -13.19784 3.880448 -3.40 0.001 -20.80646 -5.589225
3 | -11.14087 4.189539 -2.66 0.008 -19.35553 -2.926199
|
age | -.4751707 .0630779 -7.53 0.000 -.5988511 -.3514903
|
group#c.age |
2 | .2484541 .0887425 2.80 0.005 .0744516 .4224565
3 | .2903105 .1107395 2.62 0.009 .0731774 .5074437
|
_cons | 90.55752 3.009782 30.09 0.000 84.65606 96.45897
------------------------------------------------------------------------------
. margins group, dydx(age) pwcompare(pv)
Pairwise comparisons of average marginal effects
Model VCE : OLS Number of obs = 3,000
Expression : Linear prediction, predict()
dy/dx w.r.t. : age
-----------------------------------------------------
| Contrast Delta-method Unadjusted
| dy/dx Std. Err. t P>|t|
-------------+---------------------------------------
age |
group |
2 vs 1 | .2484541 .0887425 2.80 0.005
3 vs 1 | .2903105 .1107395 2.62 0.009
3 vs 2 | .0418565 .1103668 0.38 0.705
-----------------------------------------------------
R (emmeans)
Call:
lm(formula = y ~ group * age, data = m)
Residuals:
Min 1Q Median 3Q Max
-75.335 -14.841 -0.036 14.469 73.169
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 90.55752 3.00978 30.088 < 2e-16 ***
group2 -13.19784 3.88045 -3.401 0.00068 ***
group3 -11.14087 4.18954 -2.659 0.00787 **
age -0.47517 0.06308 -7.533 6.52e-14 ***
group2:age 0.24845 0.08874 2.800 0.00515 **
group3:age 0.29031 0.11074 2.622 0.00880 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21.19 on 2994 degrees of freedom
Multiple R-squared: 0.03412, Adjusted R-squared: 0.03251
F-statistic: 21.15 on 5 and 2994 DF, p-value: < 2.2e-16 contrast estimate SE df t.ratio p.value
1 - 2 -0.2485 0.0887 2994 -2.800 0.0051
1 - 3 -0.2903 0.1107 2994 -2.622 0.0088
2 - 3 -0.0419 0.1104 2994 -0.379 0.7045 Plotting
Code
Stata
. regress y i.group##c.age
Source | SS df MS Number of obs = 3,000
-------------+---------------------------------- F(5, 2994) = 21.15
Model | 47475.1592 5 9495.03183 Prob > F = 0.0000
Residual | 1343957.85 2,994 448.883716 R-squared = 0.0341
-------------+---------------------------------- Adj R-squared = 0.0325
Total | 1391433.01 2,999 463.965657 Root MSE = 21.187
------------------------------------------------------------------------------
y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
group |
2 | -13.19784 3.880448 -3.40 0.001 -20.80646 -5.589225
3 | -11.14087 4.189539 -2.66 0.008 -19.35553 -2.926199
|
age | -.4751707 .0630779 -7.53 0.000 -.5988511 -.3514903
|
group#c.age |
2 | .2484541 .0887425 2.80 0.005 .0744516 .4224565
3 | .2903105 .1107395 2.62 0.009 .0731774 .5074437
|
_cons | 90.55752 3.009782 30.09 0.000 84.65606 96.45897
------------------------------------------------------------------------------
. margins group, at(age = (20(10)60))
Adjusted predictions Number of obs = 3,000
Model VCE : OLS
Expression : Linear prediction, predict()
1._at : age = 20
2._at : age = 30
3._at : age = 40
4._at : age = 50
5._at : age = 60
------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_at#group |
1 1 | 81.0541 1.793005 45.21 0.000 77.53845 84.56975
1 2 | 72.82534 1.284642 56.69 0.000 70.30647 75.34421
1 3 | 75.71945 1.27105 59.57 0.000 73.22723 78.21167
2 1 | 76.30239 1.219243 62.58 0.000 73.91176 78.69303
2 2 | 70.55818 .8030164 87.87 0.000 68.98366 72.1327
2 3 | 73.87085 .8136044 90.79 0.000 72.27557 75.46613
3 1 | 71.55069 .744313 96.13 0.000 70.09127 73.0101
3 2 | 68.29101 .6470303 105.55 0.000 67.02234 69.55968
3 3 | 72.02224 1.168425 61.64 0.000 69.73125 74.31324
4 1 | 66.79898 .6459221 103.42 0.000 65.53249 68.06548
4 2 | 66.02384 .9857706 66.98 0.000 64.09099 67.9567
4 3 | 70.17364 1.93012 36.36 0.000 66.38915 73.95814
5 1 | 62.04727 1.037397 59.81 0.000 60.01319 64.08136
5 2 | 63.75668 1.517933 42.00 0.000 60.78038 66.73298
5 3 | 68.32504 2.782515 24.56 0.000 62.86921 73.78088
------------------------------------------------------------------------------
. marginsplot, recastci(rarea) ciopts(fcolor(%25) lcolor(%25))
Variables that uniquely identify margins: age group
The recastci(rarea) option plots the confidence region, rather than just confident intervals at each estimated point (to more closely mimic R’s output). The ciopts() allows you to pass options for the confidence intervals; these particular options just make the confidence region (and it’s outline) transparent.
R (interactions)
Call:
lm(formula = y ~ group * age, data = m)
Residuals:
Min 1Q Median 3Q Max
-75.335 -14.841 -0.036 14.469 73.169
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 90.55752 3.00978 30.088 < 2e-16 ***
group2 -13.19784 3.88045 -3.401 0.00068 ***
group3 -11.14087 4.18954 -2.659 0.00787 **
age -0.47517 0.06308 -7.533 6.52e-14 ***
group2:age 0.24845 0.08874 2.800 0.00515 **
group3:age 0.29031 0.11074 2.622 0.00880 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21.19 on 2994 degrees of freedom
Multiple R-squared: 0.03412, Adjusted R-squared: 0.03251
F-statistic: 21.15 on 5 and 2994 DF, p-value: < 2.2e-16